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1.  Introduction 


In  multi-physics  codes  simulating  large  material  flows,  representations  of  the  un¬ 
derlying  material  structure  must  transport  through  the  computational  grid.  Material 
enters  and  leaves  computations  cells.  Quantities  like  energy,  density,  and  material 
state  are  determined  by  a  weighted  average  of  the  state  entering  a  cell  and  what 
remains  after  material  exits.  Errors  are  introduced  by  ad  veering  structure  variables 
through  the  grid,  and  algorithms  minimizing  these  errors  improve  the  quality  of  the 
solution.  For  solid  mechanics,  one  of  the  more  challenging  quantities  to  advect  is  the 
deformation  gradient  tensor  because  nonlinear  combinations  of  the  9  components 
correspond  to  physical  quantities.  For  example,  the  determinant  of  the  deformation 
gradient  describes  the  volume  change.  There  is  little  chance  that  a  linear  combina¬ 
tion  of  the  9  individual  components  of  2  deformation  gradients  will  yield  the  same 
determinant  as  taking  the  linear  combination  of  the  determinants. 

One  strategy  to  reduce  advection  error  is  to  decompose  the  deformation  gradient, 
F,  into  smaller,  physically  meaningful  pieces  and  to  advect  those.  The  multi-step 
decomposition  being  investigated  at  the  US  Army  Research  Faboratory  (ARE)  be¬ 
gins  with  an  F  =  RU  decomposition  of  the  deformation  gradient,  where  R  is  an 
orthogonal  rotation  matrix  and  U  is  the  symmetric  right  stretch  tensor.  The  focus 
of  this  technical  brief  is  on  advection  of  the  rotation. 

Rotation  of  an  object  or  a  micro  structure  can  be  represented  in  many  forms,  such  as 
Euler  angles,1  angle-axis  pairs,  Rodrigues  vectors,2  quaternions  and  rotation  matri¬ 
ces.3’4  While  3  parameters  are  sufficient  to  define  a  rotation,  as  in  Euler  angles,  more 
parameters  are  often  used  for  convenience  or  other  reasons.  Consider,  for  example, 
an  angle-axis  pair,  which  is  defined  in  terms  of  a  rotation  u  about  the  unimodular 
axis  c.  While  there  appear  to  be  4  parameters,  the  axis  is  a  unit  vector,  so  there  is  an 
additional  constraint,  leaving  only  3  independent  parameters.  Likewise,  a  quater¬ 
nion  representation  has  4  parameters  and  a  constraint  equation.  A  rotation  matrix 
has  9  parameters,  with  6  constraints  from  the  orthogonality  requirement.  Advect- 
ing  large  numbers  of  redundant  parameters  is  not  ideal  since  they  must  be  adjusted 
consistently  after  advection  to  satisfy  the  constraint  equations. 

From  the  perspective  of  having  no  redundant  parameters,  Euler  angles  and  Ro¬ 
drigues  vectors  both  appear  suitable.  However,  Euler  angles  are  not  periodic  with 
continued  rotation.  A  body  is  in  an  identical  orientation  after  a  360°  rotation  about 
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an  axis,  but  the  Euler  angles  do  not  naturally  return  to  the  same  values.  A  step 
change  in  angles  is  necessary  to  keep  the  angles  within  a  fundamental  range.  The 
lack  of  periodicity  is  not  an  issue  for  many  applications,  but  advection  requires  a 
continuous  parameterization,  so  Euler  angles  are  not  suitable  for  current  consid¬ 
erations.  Rodrigues  vectors,  on  the  other  hand,  use  periodic  functions.  However, 
the  components  are  unbounded  as  the  cosine  function  in  the  denominator  passes 
through  zero.  A  singularity  cannot  be  accommodated  so  Rodrigues  vectors  are  also 
not  suitable  in  the  current  context. 

Angle-axis  pairs  and  quaternions  both  have  one  redundant  parameter.  Similar  to 
Euler  angles,  the  angle-axis  pairs  suffer  from  inability  to  represent  periodicity  in 
rotations.  The  angle  continuously  increases  with  rotation.  Quaternions  are  built 
on  periodic  functions  and  are  bounded.  In  addition,  errors  introduced  by  noise  in 
quaternions  are  less  severe  than  for  other  methods,5  such  as  Euler  angles,  so  there 
is  additional  incentive  for  their  use.  Methods  based  on  quaternions  will  be  explored 
in  further  detail  in  this  report. 

2.  Quaternion  Representation  of  a  Rotation  Matrix 

The  rotation  of  a  body  can  be  represented  in  terms  of  a  rotation,  to,  about  a  uni- 
modular  axis,  c.  Quaternions  are  related  to  angle-axis  pairs,  but  they  employ  cyclic 
functions  that  are  able  to  capture  periodicity.  In  terms  of  the  rotation  angle,  u,  and 
axis,  c,  the  4  components  of  the  quaternion  are  defined  in  a  rectangular  Cartesian 
coordinate  system  as 

A  =  ci  sin(a;/2) 

H  =  c2  sin(o;/2) 

...  J  v-U 

v  —  c3sin(cj/2) 
p  =  cos(o;/2) 

subject  to  the  constraint 

A2  +  p2  +  z/2  +  p2  =  1  .  (2) 
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Quaternions  can  be  used  to  construct  the  rotation  matrix,  R,  uniquely  by 


Rij 


A2  —  p2  —  u2  +  p 2  2(A/i  —  up)  2(u\  +  pp) 

2(Ap  +  i/p)  p2  —  i/2  —  A2  +  p 2  2  (pi/  —  Ap) 

2(i/ A  —  pp)  2  (pi/  +  Ap)  i/2  —  A2  —  p2  +  p2 


(3) 


When  determining  quaternions  from  a  rotation  matrix,  there  is  a  nonuniqueness  as 
2  sets  of  quaternions  can  be  determined  for  any  general  rotation  matrix.  Inspection 
of  Eq.  3  reveals  that  a  sign  indeterminacy  stems  from  the  quaternions  only  appear¬ 
ing  as  products.  Changing  the  sign  on  all  of  the  quaterions  gives  the  same  rotation. 
Redundancy  is  also  expected  because  of  the  half  angle  in  Eq.  1;  it  takes  2  full  rota¬ 
tions  to  run  the  range  of  the  periodic  functions.  A  common  means  of  dealing  with 
this  problem  is  to  ensure  that  one  of  the  quaterion  components  is  always  positive. 
While  this  guarantees  a  unique  solution  acceptable  for  many  applications,  it  forces 
a  discontinuous  sign  change  in  other  components.  The  discontinuity  could  create 
issues  if  quaternions  are  used  to  advect  a  rotation. 


Figure  1  shows  examples  of  quaternion  components  for  400°  rotations  of  2  different 
initial  rotation  matrices.  The  cosine  component,  p,  was  forced  to  be  positive.  As  p 
hits  zero,  where  there  would  have  been  a  sign  change,  other  quaternion  components 
flip  sign  creating  large  discontinuities. 


Fig.  1  Evolution  of  quaternion  components  for  2  different  initial  rotation  matrices 

Had  another  of  the  quaternion  components  been  chosen  to  be  always  positive,  the 
discontinuity  would  have  occurred  in  another  place. 
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3.  Continuity  from  Evolution 


Continuity  conditions  can  be  applied  to  evolve  the  quaternions  if  the  rotations  are 
associated  with  the  deformation  gradient.  The  deformation  gradient  begins  as  the 
identity  and  evolves  continuously  over  time.  Generally,  the  deformation  gradient  at 
the  beginning  of  the  time  step  will  be  available  when  the  deformation  gradient  at 
the  end  of  the  step  is  calculated.  This  allows  selection  of  the  sign  on  the  quaternion 
to  keep  the  values  continuous.  Within  a  finite  element  code,  as  long  as  the  deforma¬ 
tion  gradients  in  neighboring  elements  continue  to  evolve  smoothly,  discontinuities 
should  never  exist  between  adjacent  elements.  Barton6  has  used  the  sign  of  the  inner 
product  of  the  beginning-of-step  and  end-of-step  quaternions  to  trigger  a  sign  flip. 

Pop  +  A0A  +  /Jo/J  +  u0u  <  0  .  (4) 


If  the  inner  product  is  negative,  the  sign  of  all  of  the  quaternion  components  is 
changed.  The  rotations  shown  in  Fig.  1  are  recalculated  using  sign  continuity,  and 
the  results  are  shown  in  Fig.  2.  The  quaternions  are  now  continuous  throughout  the 
range. 


Fig.  2  Evolution  of  quaternion  components  for  2  different  initial  rotation  matrices  using  con¬ 
tinuity  from  previous  value 


This  continuity  criterion  is  particularly  well  suited  for  the  deformation  gradient, 
where,  initially,  p  =  1  and  A  =  p  =  v  =  0.  The  sum  in  Eq.  4  is  dominated  by  p 
in  the  vicinity  of  the  initial  orientation,  and  no  sign  change  will  be  indicated.  The 
remaining  values  transition  smoothly  from  positive  to  negative.  Other  quaternion 
components  dominate  Eq.  4  at  large  rotations  to  maintain  continuous  values. 
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While  the  method  ensures  continuity,  not  all  symmetries  are  captured.  Consider  2 
material  points  with  common  material  directors.  After  rotating  each  180°  in  op¬ 
posite  directions,  they  will  again  have  a  common  orientation.  However,  the  sine 
portion  of  the  quaternions  in  Eq.  1  will  have  opposite  algebraic  signs  and  will  not 
show  the  proper  orientation  relationship.  The  half  angles  create  a  parameter  space 
that  is  twice  the  size  needed,  and  the  functions  do  not  cycle  back  on  themselves  at 
the  correct  rate. 

The  criterion  also  does  not  resolve  ambiguities  in  other  initial  rotations.  Say  the 
initial  rotation  is  180°  about  the  z-axis.  The  first  2  diagonal  elements  of  the  rotation 
matrix  are  -1  and  the  off-diagonals  are  zero.  Perturbations  on  one  side  of  this  have 
the  product  pv  being  positive  and  the  product  is  negative  for  perturbations  on  the 
other  side.  The  signs  of  p  and  v  individually  are  not  determined  without  introduc¬ 
tion  of  additional  rules.  If  p  is  assumed  positive,  as  in  the  case  of  the  initial  rotation 
being  the  identity,  the  v  will  be  discontinuous  between  positive  and  negative  rotat¬ 
ing  points.  Alternative  rules  would  have  to  be  applied  in  such  situations. 

It  is  not  anticipated  that  the  deficiencies  noted  above  will  be  an  issue  for  the  defor¬ 
mation  gradient,  but  such  situations  may  arise  for  the  lattice  orientations  in  crystal 
plasticity  simulations.  Alternative  methods  may  be  useful  for  the  more  general  case. 

4.  Proposal  for  a  Unique  Representation 

The  solution  proposed  here  is  to  retain  the  quaternion  products  from  the  rotation 
matrix  so  that  unique  determination  of  the  algebraic  sign  is  not  necessary.  Hence, 
the  products  A p,  pp,  vp,  and  p2  are  advected.  There  is  an  obvious  issue  when  p  —  0 
and  all  of  the  component  products  to  be  advected  are  zero.  The  proposed  solution 
includes  a  lower  limit  on  p  such  that  p2  >  e2  for  some  small  e.  With  this  condition, 
the  quantities  advected  are  A p,  pp,  vp,  and  pp,  where 


p  if  p2  >  e2 


(5) 


P 


e  otherwise 


Given  a  rotation  matrix,  the  procedure  is  to  first  determine 


(6) 
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and  compare  the  value  with  e2.  If  p2  >  e2,  then 

1 

Ap  =  -  (-R32  —  R23) 
PP  —  ^  (R13  —  R31) 
vp  =  |  (-R21  -  ^12) 
PP  =  P2 


(7) 


Ifp2  <  e2  the  rotation  matrix  diagonal  is  used  to  compute  components  as: 

Ap  =  -e  sign  (f?32  —  -R23)  a/  +-Rh  —  R22  ~  R33  +  1 

pp=\e  sign  (f?i3  -  R31)  \/ —R\i  +  R22  ~  R33  +  1 

2  (8) 

z/p  =  -e  sign  (f?21  -  f?i2)  A/-i?n  -  i?22  +  #33  +  1 
PP  =  e 

This  latter  case  only  applies  when  rotation  angles  are  near  180°,  where  A2  +  p2  + 
v2  ~  1  and  the  magnitude  of  the  advected  quantities  is  set  by  the  parameter  e. 

Following  advection,  limits  should  be  imposed  on  the  advected  values  and  the  con¬ 
straint  equation  must  be  applied  to  ensure  valid  values  and  a  pure  rotation.  Since 
each  of  the  quaternion  factors  is  bounded  between  -1.0  and  1 .0,  the  advected  quater¬ 
nion  products  (Ap,  pp,  up,  and  pp)  are  also  bounded  between  -1.0  and  1.0.  Appli¬ 
cation  of  the  constraint  from  Eq.  2  determines  the  scaling  factor,  p*2,  through 

p*2  =  (Ap)2  +  (pp)2  +  (up)2  +  (pp)2  .  (9) 


Dividing  through  Eq.  9  by  p*2,  it  can  be  seen  that  the  ratio  p2/p*2  multiplies  each 
of  the  quaternion  values  to  satisfy  the  constraint.  If  there  are  no  advection  errors  the 
ratio  is  1.0.  The  quaternion  products  needed  to  reconstruct  the  rotation  matrix  can 
then  be  obtained  without  loss  of  sign  information  from  the  advected  quantities  and 
the  scaling  factor,  p*2,  using 


A2  = 


(Ap)' 

p*2 


Ap  = 


Ap  pp 


P 


and 


P2  = 


p*2 


(10) 
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Similar  expressions  are  applied  to  determine  the  remaining  terms.  These  quantities 
can  be  used  directly  in  Eq.  3  to  calculate  the  rotation  matrix.  The  value  of  p*2  will  be 
nonzero  unless  advection  error  results  in  all  4  advected  quaternion  products  being 
zero.  The  limiting  value,  e,  was  introduced  to  prevent  a  zero  solution,  but  it  should 
be  guarded  against  in  the  algorithm  nonetheless.  This  can  be  ensured  by  enforcing 
p*2  >  e2  when  evaluating  Eq.  9  by  scaling  each  of  the  terms  on  the  right  hand  side. 
It  is  critical  that  the  values  used  in  Eq.  10  satisfy  Eq.  9. 


Using  the  representation  in  Eq.  7  and  Eq.  8,  the  4  advected  components  are  plotted 
in  Fig.  3  for  the  same  rotations  as  in  Fig.  1.  It  can  be  seen  that  the  functions  are 
smooth  and  continuous.  One  of  the  points  in  each  figure  has  pp  =  0  identically,  and 
the  other  quantities  are  nonzero  because  of  the  e  offset. 


Fig.  3  Evolution  of  quaternion  products  for  2  different  initial  rotation  matrices 

It  is  interesting  that  the  quaternion  products  proposed  create  a  representation  in 
terms  of  the  full  rotation  angle  from  an  angle-axis  pair  rather  than  the  half  angle 
shown  in  Eq.  1.  The  products  are: 

A p  =  Ci  sin(o;/2)  cos(cu/2) 
pp  =  C2sin(o;/2)  cos(cu/2) 
up  =  c3  sin(o;/2)  cos(cu/2) 
pp  =  cos2(cu/2)  =  ^  [1  +  c 


=  -  ci  sm(u;J 
=  ^  c2  sin  (a;) 
=  \°3  sin(w) 


(U) 


3os(o;)] 
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Here  the  angular  range  for  periodicity  corresponds  to  the  physical  rotation.  The 
angular  range  for  periodicity  for  quaternions,  on  the  other  hand,  is  twice  that  of 
the  physical  rotation.  Keeping  the  correct  angular  range  eliminates  the  antipodal 
problem.  However,  the  representation  in  Eq.  1 1  does  not  have  a  simple  form  for  the 
constraint  equation,  as  Eq.  9  involves  both  p 2  and  p4 . 

The  troublesome  special  cases  described  with  the  continuity  approach  in  Section  3 
do  not  create  the  same  issues  with  the  proposed  representation.  Since  the  full  ro¬ 
tation  angle  is  used  in  Eq.  1 1  rather  than  the  half  angle,  2  material  points  rotated 
180°  in  opposite  directions  will  have  the  same  representation.  In  addition,  initial 
orientations  of  180°  are  zero,  so  there  will  not  be  discontinuous  behavior  at  these 
points. 

The  potential  issue  with  the  proposed  representation  is  the  method  to  avoid  the 
singularity  when  p  is  zero.  This  will  require  further  exploration  and  use  to  uncover 
any  difficulties. 

5.  Summary  and  Conclusions 

Two  methods  are  explored  for  dealing  with  the  antipodal  ambiguity  in  represent¬ 
ing  rotations  by  quaternions.  Enforcing  continuity  between  beginning  and  end  step 
values  should  work  well  for  rotations  associated  with  the  deformation  gradient.  An 
alternative  method  based  on  the  sine  of  the  full  rotation  angle  is  proposed  for  situa¬ 
tions  where  continuity  cannot  be  enforced  or  for  more  arbitrary  initial  orientations. 
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